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The  power  series  method  developed  by  the  author  II]  is  applied  to  the  heat 
equation.  Highly  accurate  semi-discrete  systems  of  equations  in  t and  in  x 
are  generated  and  are  made  stable  by  proper  choice  of  parameters.  A totally  dis- 
crete scheme  is  produced  that  represents  arbitrarily  high  accuracy  in  both  x and 
t.  Stability  analysis  indicates  that  while  arbitrary  order  in  t may  be  stable, 
the  order  of  accuracy  in  x is  restricted  to  be  less  than  16  and  certain  geometri- 
cal restrictions  on  the  step  sizes  must  be  met.  Truncation  errors  are  examined 
and  a consistency  condition  is  obtained  that  further  restricts  the  step  sizes. 

The  scheme  is  shown  to  coincide  with  Keller's  Box  Scheme  14]  in  its  lowest  order. 
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SIGNIFICANCE  AND  EXPLANATION 

In  MRC  TSR  # 1923  the  author  introduced  a method  for  solving  partial  differen- 
tial equations  by  substitution  of  power  series  into  the  differential  equation.  In 
that  report  some  technical  properties  of  the  method  were  discussed  in  connection 
with  ordinary  differential  equations.  This  paper  presents  the  first  application 
of  the  method  to  a partial  differential  equation  and  takes  as  a basic  example  the 
heat  equation.  The  method,  in  effect,  substitutes  a double  power  series  with  vari- 
ables X and  t,  directly  into  the  equation  and  finds  the  correct  coefficients 
to  satisfy  the  differential  equation  and  the  boundary  conditions  in  a grid  of  cells. 

This  simple  procedure  encounters  many  complications,  however,  and  the  bulk  of  the 
paper  involves  stepping  slowly  toward  the  goal  of  obtaining  equations  that  govern 
the  coefficients.  Conditions  on  the  built  in  parameters  are  derived  that  guarantee 
the  errors  due  to  the  approximation  to  be  small  when  introduced  and  not  to  grow  i 

with  time.  ’ 

I 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary 
lies  with  MRC,  and  not  with  the  author  of  this  report. 
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POWER  SERIES  METHODS  II  - THE  HEAT  EQUATION 


Robert  D.  Small 


Introduction 

In  a previous  paper  [11,  hereafter  referred  to  as  I,  this  author  proposed  a method  whereby 
power  series  are  used  to  generate  arbitrarily  accurate  finite  difference  schemes  for  initial 
value  problems  in  ordinary  differential  equations.  It  was  shown  t.here  that  a scheme  may  be 
made  A-stable  in  the  sense  of  Dalquist  [2]  by  choosing  the  point  of  expansion  of  the  power 
series  to  be  the  center  of  the  interval  in  which  the  solution  is  sought.  This  paper,  the 
sequel  to  I,  applies  the  method  to  a pcurticular  partial  differential  equation,  the  heat  equa- 
tion, to  illustrate  in  detail  how  a stable  scheme  with  very  high  accuracy  in  both  variables  of 
a partial  differential  equation  may  be  attained.  Since  we  shall  rely  heavily  upon  the  work  of 
I we  refer  to  equations  of  that  paper  with  a preceding  I and  similarly  a I precedes  the  figure 
number. 

In  a manner  analogous  to  that  of  I,  we  fit  power  series  solution  surfaces  over  the  domain 
that  has  been  partitioned  into  rectangular  cells.  First  infinite  power  series  are  found,  and 
after  initial  and  boundary  conditions  are  satisfied  the  series  are  truncated.  The  truncations 
that  lead  to  stable  schemes  are  not  as  straightforward  as  in  the  case  of  ordinary  differential 
equations.  We  find  semi-discrete  schemes,  discretized  in  t in  section  2 and  in  x in  section 
3.  When  these  schemes  are  made  stable  by  proper  choice  of  some  of  the  built  in  parameters, 
comparison  of  the  two  schemes  leads  to  the  totally  discrete  scheme  of  section  4.  The  stability 
analysis  of  this  scheme  then  relates  back  to  the  analysis  of  I.  Section  5 deals  with  truncation 
error  and  consistency  and  in  section  6 we  note  briefly  that  the  lowest  order  scheme  coincides 
with  Keller's  Box  Scheme. 
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1.  Set-up  of  the  Problem 

We  now  apply  the  power  series  method  to  the  heat  equation 

(1,  • 

3x 

The  domain  is  partitioned  into  a grid  of  cells  of  length  2h  in  the  x-direction  and  height  2k 
in  the  t-direction  as  illustrated  in  Figure  1.  Each  cell  is  identified  by  an  ordered  pair  of 
integers  (i,j)  indicating  its  position  in  the  grid.  The  solution  or  its  approximation  defined 
there  is  denoted  u^^(x^,t^).  Co-ordinates  within  each  cell  are  denoted  by  (x^,t^)  and  are 
measured  with  respect  to  an  origin  that  lies  on  the  vertical  mid-line  of  the  cell  a distance 

2Xh  from  the  bottom.  They  are  related  to  the  original  ones  by 

x^  = X - (2i-l)h 
tj  = t - 2(j-l+X)k  . 

This  choice  is  made  for  comparison  with  the  treatment  of  ordinary  differential  equations  in  I 
and  X will  later  be  taken  to  be  1/2  for  reasons  similar  to  those  in  I.  For  illustration 

we  take  a rectangular  grid  of  L cells  in  width  in  the  x-direction. 

The  solution  u is  prescribed  as  U(x)  initially,  g^^Ct)  on  the  left  side  of  the  domain 
and  gjCt)  on  the  right.  Then  if  the  domain  is  0 x £ 2Iih,  0 t < ■»,  we  can  express  all 
variables  in  local  co-ordinates  giving  initial  cmd  boundary  conditions  as 

Uii(Xi,-2Xk)  = i = 1,2,'”,L 


j = l,2,”-,“ 


where  the  quantities  on  the  right  side  are  defined  by 


^i**i'  ” U(x),  gj^^(t^)  = gj^(t),  = g2(t) 

Other  domain  shapes  and  boundary  conditions,  including  free  surfaces  Ceui  be  handled  without 
major  modification  of  this  basic  method. 


and 


For  interior  cell  boundaries  u is  continuous  across  the  horizontal  boundaries  and  u 


3u 

3x 


are  continuous  across  vertical  boundaries.  These  conditions  are 


(3) 


Uij (x^,-2Xk) 


u..(-h,t.) 


3u.  . 


u.  . , (x.  ,2(1-X)k) 


= u.  , .(h,t.) 
1-1.3  3 


3u.  . 

= "3^ 


fi  = 1,2,---,I 
= 2,3,-'-,“ 


1 — 2 f3  f * * * flj 


l,2r’ 


Commas  are  used  to  separate  subscript  expressions  for  clarity  and  are  omitted  when  no  ambiguity 
arises.  No  two  subscripts  are  multiplied  and  subscripts  are  not  used  to  indicate  partial 
differentiation. 


We 


determine  U..(x,,t.)  within  each  cell  by  finding  an  infinite  power  series  in  x.  and 
1]  1 3 1 

t^  satisfying  (1)  and  then  we  truncate  it.  The  process  of  truncating  in  such  a way  as  to 
leave  a stable  and  well  posed  scheme  is  very  complicated  and  it  is  much  simplified  by  treating 
one  variable  first  and  then  beginning  again  with  the  other.  We  encounter  conditions  that  must 
be  satisfied  for  each  semi-discrete  scheme  to  be  stable  and  then  we  combine  the  two  sets  of 
conditions  in  the  totally  discrete  formulation. 


2.  Discretization  in  t 


Similarly,  substitution  of  (4)  into  (3)  gives 


I " - (-2Xk)"  = y 

£.  n!  ^ 

n=0  n= 


J2n)  , ^ 

" a . . (x . ) 

y — — (2(i-x)k)" 

n=0  "• 


i = 1,2,”',L 


I \ I \ I 1=2, 3, •'•,L 

a.,  (-h)  = a._  (h)  < ] = 2, 3, •••,“> 

^ I n = 0,1 , • • • 


We  see  now  that,  for  fixed  j and  assuminq  that  the  a.  . are  known  functions,  (7)  is 

1 , 3“1 

a coupled  _jt  of  L infinite  order  ordinary  differential  equations  in  the  L functions  a^^ 
and  there  are  the  equivalent  of  L sets  of  boundary  data  given  by  (6)  and  (8) . In  addition, 
(5)  is  a set  of  ordinary  differential  eouations  that  starts  the  whole  process  at  j = 1.  This 
is  an  unnecessarily  complicated  interpretation,  however,  since  (8)  states  that  a^^  is  con- 
tinuous to  all  orders  of  derivatives  in  the  x-direction  and  hence  we  may  disregard  the  sub- 
script i and  solve  for  one  function  a^  in  a horizontal  strip  of  cells  at  time  step  j that 
extends  the  entire  width  of  the  domain. 

So  far  the  problem  is  exact.  We  may  approximate  it  by  truncating  the  series  in  (5)  and 
(7)  at  n = N which  leaves  an  order  2N  ordinary  differential  equation  for  the  one  function 
a^^y  Then  the  proper  boundary  conditions  selected  from  (6)  are  those  that  apply  for 
n = 0,1,"**,N-1,  and  the  continuity  conditions  (8)  apply  for  the  range  n = 0,1, • • • ,2N-1. 
Having  solved  for  all  the  a^^  one  obtains  the  solution  by  substituting  into  (4)  which 

has  been  truncated  at  n = N.  Thus  the  semi-discrete  scheme  of  order  N,  discretized  in  t. 


N a<f>(x.) 
r il  1 

n^0~^ 


— (-2Xk)"  = U.(x.)  i = 1,2,--”,I 


"ir  <-h) = (0) 


(2n)  ...  (n) 

\j  = ^2j 


j = 1,2,- 


n = 0,1,*”,N-1 


H a‘^"’(x.l  N a!^"’  (x.) 

J ,-2Ak)"  = y ^ 

n=0  n=0  ■ 


(2(1-A)k)‘ 


i — 1,2,*"“,L 


j = 2,3,' 


i = 2, 3, • • • ,L 


(n)  (n)  ‘ , 

a . (-h)  = a.  , , (h)  < i = 2,3,' 

11  1-1-1  \ 


n = 0, 1, • • • ,2N-1 


X.  ^2n)  , ^ 

N a..  (x.) 


u..(x.,t.)  = I t" 


i = 1,2, •••,L 


j = 1,2,- 


The  problem  is  now  discrete  in  the  time  direction,  j having  replaced  t^ , but  still  con- 
tinuous in  space  as  x^  is  present  and  i is  superfluous.  This  system  may  be  tested  for 
stability  using  Von  Neumann's  criterion.  Since  <9)  is  a linear  system,  on  an  infinite  domain 
the  error  satisfies  the  same  equations  as  the  solution  and  we  can  examine  one  Fourier  component 
of  the  error  by  setting 

• • 

a,  . (x.)  = 6^e  ^ a 

13  1 

where  a is  any  real  number  corresponding  to  the  wave  number  of  an  error,  a is  any  amplitude 
of  error,  B gives  the  rate  of  growth  of  the  error,  which  is  to  be  determined,  and  1 = 
Substituting  this  into  (9c) , we  obtain,  for  a ^ 0, 

(10)  I ,2a\)"  = 0 . 

n=0 


Solving  for  6 we  have 


N T / / N / 

= I I2(X-l)a^)c]'/n!  / I (2Aa^k)7n!  . 

n=0  / n=0 


For  a real,  8 is  real  and  for  the  scheme  to  be  stable  we  must  have  |8l  ^ 1 so  that  the 
error  does  not  grow.  We  shall  have  stability,  therefore,  if  for  |6|  > 1,  (10)  has  no  real 
roots  for  a.  This  is  the  case  if  the  term  in  square  brac)cets  has  consistent  sign  for  all  n, 
•■M.ich  occurs  when  X > 1/2.  Thus  the  scheme  (9)  is  stable  if  we  choose  X ^ 1/2. 


3.  Discretization  in  x 


We  now  discretize  in  the  x-direction  to  obtain  further  stability  conditions  that  must  be 


met  in  the  full  discretization.  The  power  series  in  x^  that  satisfies  (1)  is 


— (m)  _ , 
“ a;  . (t .) 


(ID 


„ (x  t ) = y -ij x^”  + y -ij L „2m+i 

ij‘  i'  j'  (2m)!  i (2m+l)  ! i 

m— 0 m— 


m=0 


where  the  a.  . and  the  b.  . are  functions  of  t.  yet  to  be  determined.  Substitution  of  (11) 


ID 


into  (2)  and  (3)  gives 


a!""'  (-2Xk)  = u!^""  (0) 
ll  1 


(12a) 


b!^^  (-2X)c)  = (0) 


(12b) 


aiT’ (t J 


y . y 

^ (2m) ! " i 


bi”>(tj 


, 2m+l 


m=0 


m=0 


(2m+l)  ! 


= g,.(t.) 


liT’  (tj 


> j = 1,2,' 


y -iJ i_  + y 3-  -2m+i 

^ (2m) ! (2m+l) ! 


n¥=0 


itpsQ 


(12c) 


a!'"N-2X]c)  = a!'"'.  (2(1-X)!!;) 

1 3 1 / j-1 


bf'!'’(-2X)c)  = b!”*^  ,(2(1-X)]c) 
t3  t»j“l 


m = 0,1, • • • ,“ 
i = 1,2, • • • ,L 
j = 2,3,  • • • ,"» 


» af’!’*  (t.) 

y ^-h^" 

mio 


^‘""(t.) 


_ y L 

(2m+l)  ! 


m-O 
— (m) 


(tj 


= y i-i.j'  3 

^ (2m)! 


2Tn 


m=0 


r"(n\) 

« b (t  ) 

y i,2mtl 


m=0 


(2m+l)  ! 


(12d)< 


-I 


b!T>  (tj 


m=0 


(2m+l)  ! 


y _ij L-h^" 

^ (2m)  ! 


m=0 
— (m+1) 


a ' ( t ) 

= y ,2m4l 

^ (2m+l) ! 


” bf  .(t.) 

y 3-J-.>.3 3_  h^m 

(2m)! 


i = 2,3,- 
j = 1,2,- 


’ • ,L 


We  now  interpret  (12b)  and  (12d)  as  2L  coupled  infinite  order  ordinary  differential  equations 

in  the  2L  functions  a. . ^md  b. . with  2L  sets  of  data  in  (12a)  for  the  case  j = 1 and 

ID  13 

2L  sets  of  data  in  (12c)  for  j i 2.  The  continuity  conditions  (12c)  , however,  indicate  that 

both  a . and  b. . are  continuous  to  all  orders  of  derivatives  in  the  t-direction  and  hence 
13  ID 

j plays  no  role  in  the  solutions.  We  are  really  solving  for  2L  functions  in  L vertical 
strips  that  extend  infinitely  forward  in  time. 

We  approximate  this  exact  system  by  truncating  the  series  in  (12d)  at  finite  values  of  m. 
Displaying  only  the  summation  signs,  we  modify  (12d)  to 

«1  ”2  «1  «2 

I - I = I + I 

m=0  m=0  m=0  m=0 

(13) 

«3  ”4  «3  »4 

- I + I = I - I 

m=0  m=0  m=0  m=0 

V 

Then  both  of  equations  (12b)  are  truncated  as  in  the  first  of  (13)  so  that  equations  that  main- 
tain continuity  of  the  function  value  are  of  consistent  accuracy. 

Now  the  differential  equations  are  discrete  in  space  and  continuous  in  time.  We  postpone 
selecting  the  correct  boundary  conditions  from  among  (12a)  cmd  (12c)  until  the  upper  limits  of 
(13)  are  chosen  to  make  that  system  stable.  We  examine  one  Fourier  component  of  the  error  by 


setting 


6t.  ^ 

— . 3 2iiah  - 

a.  . (t.)  = e -"e  a 

13  3 


6t.  .■r.  ^ 

r » 3 2iiah  * 

b.  . (t . ) = e e b 

13  3 


for  real  arbitrary  a,  a,  b.  Substituting  (14)  into  (13)  we  obtain  two  linear  homogeneous 
equations  in  a and  b.  The  condition  for  non-trivial  solutions  is 


Ml  m M.  , m 7 ^2  2 M-  - m+1/ 

MS)  y (6h  ) r (Bh  ) -2icih  ^ r (Bh  ) r ( Bh  ) 

' (2m)!  (2m)!  ' ' (2m+l)  ! (2m+l)  ! 

in=0  m=0  111=0  m=0 


,,  -2iah, 
(1+e  ) 


The  expression 


simolifies  to  - cot  ah. 


r' 


since  a is  arbitrary,  this  expression  can  take  any  negative  value  and  we  set  it  equal  to  -c  . 

We  can  take  square  roots  of  (15)  if  we  choose  M,  = M,  and  M_  = M, . Then  we  must  have 

14  2 3 

Ml  = or  + 1 for  maximum  accuracy  since  this  does  not  retain  in  (13)  any  term  more 

accurate  than  one  omitted.  Talking  the  square  root  of  (15)  and  setting 


(16) 

and 

(17) 

we  obtain 
(18) 


2 

gh  ' z 


= lM/2]  , M2  = l(M-l)/2] 


1 2m  2 2m+l 

V z _ 5^  ^ _ 


m=0 


(2m)  ! 


m=0 


(2m+l)  ! 


which  is  precisely  equation  (1-16)  derived  in  I in  connection  with  A-stability  of  this  method 
and  M plays  the  role  of  N of  (1-16) . 

For  stability  we  require  that  Re  $ £ 0 in  order  that  errors  (14)  do  not  grow  with  time. 
Accordingly,  (16)  indicates  that  we  have  staibility  if  the  roots  of  (18)  satisfy 


(19) 


IT  3it  511  7it 

T < arg  z < -r-  or  — < arg  z < — 
4 — “"4  4 — — 4 


The  roots  of  equation  (18)  form  closed  curves  parametrized  by  c.  These  curves  originally 
appeared  in  Figure  I-l  aind  are  illustrated  later  in  Figure  2b.  Condition  (19)  is  satisfied  by 
these  curves  for  M ^ 15  and  hence  M = 15  is  the  greatest  order  of  accuracy  that  may  be  used 
stedily . 

We  now  choose  the  appropriate  boundary  conditions  from  (12a)  and  (12c)  to  suit  differential 
equations  (12b)  and  (.12d)  . The  serai-discrete  scherae  of  order  M,  discretized  in  x,  is 


(20a) 


If;”  (-21k)  = o!2n»(0, 

il  1 


0, !,•••, M2  I 

i = 1,2,-*",I 


b|™’ (-21k)  = (0)  m = 0,l,”‘,Mj^-l 


-10- 


obtaining  the  set  of  equations 


Equations  (21)  express  exactly  the  semi-discrete  scheme  continuous  in  x^  and  equations 

(22)  express  exactly  the  semi-discrete  scheme  continuous  in  t^ . The  former  is  stable  for  all 

N if  X ^ 1/2  and  the  latter  is  stable  for  M £ 15,  where  and  are  related  to  M 

by  (17).  These  become  fully  discrete  schemes  if  the  infinite  sums  are  truncated  but  it  is 
difficult  to  see  how  to  do  this  stably  on  either  scheme  in  isolation.  Each  can  serve  to  show 
how  to  truncate  the  expansion  in  one  of  either  x^  or  t^  and  it  remains  to  select  common 
variables  for  them  so  that  iDOth  truncations  can  be  applied  simultaneously. 

Comparing  (21e)  and  (22e)  we  see  that  the  expansion  variables  are  related  by 

(23)  c. . = d.  , d. . = d. . _ 

i],n  i],2n  i^n  i],2n+l 

Placing  these  into  (22)  we  obtain  equations  resembling  (21)  but  with  sums  infinite  in  n cuid 
finite  in  m.  From  these  two  sets  of  equations  we  ta)ce  all  the  finite  limits  and  obtain  the 
totally  discrete  scheme  which  follows. 


(24a) 


N d /m  = 0,1,***,M-1 

I (-2X30"  = uf"'>(0)  { 


n=0 


W- 

m 

n 


L«- 


The  last  formula  is  used  only  for  displaying  the  solution  after  the  been  confuted. 

The  last  term  in  (24e)  appears  because  d. . is  not  computed  and  must  be  deleted  from  the 

13  »M+zN 

double  sum. 

We  now  have  a finite  difference  scheme  of  order  (M,N)  for  the  heat  equation  (1) . This 
scheme  is  now  tested  for  stability  by  taking  a Fourier  component  of  the  error. 


^ i - 

to 

» 

tj?. 

I 


d.  . = 

i3n  n 


Substitution  into  (24c)  and  (24d)  gives  the  homogeneous  set  of  equations 

N / > V rt"!  < 1 \ n 

I (2k.)"d_^.,_  = 0 m = 0,1, •••,«-! 


n=0 


n! 


m+2n 


(25) 


M , ,,n  -2iah 

h^d  ^ =0  n = 0, !,•••, 2N-1  . 

m!  m+n 

m=0 


The  condition  for  non-trivial  solutions  is  that  the  following  matrix  have  zero  determinant. 


D = 


'0 

0 


0 


0 Y2  ° \ 


0 Yj^  0 


0 

6, 


Y^  ....  0 


0 Y^  0 


Vl  ° 


'N 

0 


t 


(-1)"-B"^(1-X)"  .n 

where  Y = : (2k) 

n n ! 


, ,.m  -2icih 

, . (-1)  -e  , m 

and  6 = : h 

m m! 


Matrix  D is  of  the  form  of  a Sylvester  matrix  and  its  determinant  is  called  the  result- 
ant of  the  two  polynomials 


5- 


N 2 M 

n=0  m=0 


It  may  be  shown  that  det  D = 0 if  and  only  if  Pj^(z)  and  P2(z)  have  at  least  one  common 
zero.  The  classical  proof  (see  Bocher,  [3],  for  example)  begins  with  the  assumption  that  pj^(z) 
and  P2(z)  have  a common  factor  and  builds  from  the  fact  that  there  are  no  polynomials  f^(z) 


and  f2(z)  such  that 


fl(z)Pi(z)  + f2(z)p2(z)  » 1 . 


We  present  an  alternative  proof  here  since  the  classical  one  cannot  be  extended  to  treat  the 
more  complicated  matrix  that  we  encounter  in  section  5 in  connection  with  truncation  errors . 

We  note  first  that  det  D = 0 is  an  algebraic  equation  of  degree  M in  6 for  fixed  n 
and  hence  there  are  M solutions  for  6.  Further  we  note  that  the  substitution 


(26) 


d = z 
n 


n = 0,1,*‘-,2N+M-1 


into  (25)  gives  only  two  independent  equations/ 


(27) 


Pj^(z)  = 0 , P2(z)  = 0 , 


since  all  others  differ  from  one  of  these  only  by  a factor  of  a power  of  z.  If  equations 
(27)  hold  simult2meously  then  (25)  are  satisfied.  It  remains  to  prove  that  there  are  not  fur- 


ther solutions  alternative  to  (26) . For  fixed  o,  P2(z)  = 0 gives  M solutions  for  z.  Then 


since  Pj^(z)  =0  is  linear  in  6/  the  substitution  for  z gives  M solutions  for  S as  re- 


•'PI 


quired  and  no  alternative  solutions  exist.  Thus  there  are  non-trivial  solutions  for  the  error 


if  (27)  holds 


System  (24)  is  stable  if  (27)  holds  only  for  |b1  ^ 1 and  not  otherwise.  To  investigate 
this  condition  we  assume  that  |s|  >1  and  try  to  show  that  the  zeros  of  p, (z)  and  P2(z) 


are  disjoint  for  all  real  a.  Then  since  det  D = 0 is  an  algebraic  equation  for  6 
must  be  solutions  for  6 and  they  must  satisfy  |b|  £ 1. 


We  first  examine  the  equation  p. (z)  = 0.  It  happens  to  be  a form  of  equation  (1-12) . To 


find  the  (joundary  t>etween  the  steible  and  unstable  regions  we  set  |r|  = 1 in  (1-12)  and  this 


5 


6 


2 3 4 

Figure  2a.  Roots  of  = 0 

The  values  of  N are  Indicated.  The  regions  where  |8| 


> 1 are 


shaded. 


for  real  i(i  . The  result  is 


— i A 

we  do  by  setting  R = e 


N 

y 

(l-X)"-e^*(-X)" 

L 

n=0 

n! 

By  comparison, 

if  we  wish  to  set 

1b|  = 1 in  p^(z)  ■ 

{ 

N 

y . 

(l-M"-e^^-X)"  : 

V aO/ 

L 

n=0 

n! 

The  roots  for  qh  found  in  I correspond  to  the  roots  for  2kz  . We  found  it  necessary  to  set 
X = 1/2  in  order  for  inequality  (1-14)  to  hold  for  N ^ 3 and  this  corresponds  here  to  one  of 

the  conditions  necessary  to  have  |b1  £ 1 for  all  a when  N £ 3.  Thus  setting  X = 1/2,  we 

2 

have  in  Figure  I-l  a plot  of  the  roots  of  (28)  where  l^z  is  the  variable.  For  comparison 
with  the  other  polynomial,  this  plot  is  mapped  into  the  plane  with  variable  z and  its  first 


quadrant  is  shown  in  Figure  2a.  The  regions  where  |B|  >1  appear  shaded. 
Now  considering  the  equation  = 0»  if  it  is  divided  by  e - 


1 and  the  terms 


in  a collected,  we  obtain 


1 ,2m 

y (hz)  ..  ^ I 


2 ,,  , 2m+l 

y (hz) 

(2m+l)! 


where  c = cot  ah.  We  have  seen  this  equation  before  as  (1-16)  . Its  roots  for  the  varictble 
hz,  parameterized  by  c,  form  the  closed  curves  and  the  imaginary  axis  of  Figure  I-l  and 
apply  for  orders  M rather  than  N.  The  first  quadrant  of  this  plot  is  reproduced  in  Figure 

2b. 

The  stability  of  a given  scheme  of  order  (M,N)  depends  on  the  positions  of  the  curves 
corresponding  to  N in  Figure  2a  euid  M in  Figure  2b.  If  any  curve  of  Figure  2b  intersects 
a shaded  region  of  Figure  2a  when  the  scales  are  made  consistent,  then  the  scheme  is  unstable. 
Clearly  instability  can  be  avoided  if  one  of  the  rings  of  Figure  2a  intersects  with  a ring  of 
Figure  2b  since  the  value  of  h or  k can  be  adjusted  to  move  the  rings  apart.  Instcibility 
cannot  be  avoided  when  a ring  of  Figure  2b  crosses  a 45°  line  of  Figure  2a  and  enters  the  large 
shaded  region.  This  begins  at  M = 16.  The  rings  of  Figure  2a  never  touch  the  imaginary  axis 
and  thus  arbitrarily  large  values  of  N may  be  used. 


- i 


In  summary,  the  fully  discrete  scheme  (24)  with  X = 1/2  is  unconditionally  stable  for 
l^Mf^4,  1;^N^4.  For  5 ^ M ^ 15  and  5 £ N < “ there  are  mild  restrictions  on  h and 
Ir  for  stability  that  are  found  graphically  from  Figures  2a  and  2b.  For  M 16  the  scheme  is 
unstable.  These  stability  properties  resemble  those  of  the  ordinary  differential  equation  of  I 
where,  for  N 5 the  step  size  h had  to  be  adjusted  so  that  the  values  of  ^ avoided  the 
small  regions  of  Figure  I-l. 


I 


5.  Truncation  Error  and  Consistency 


A desirable  attribute  of  a finite  difference  scheme  is  that  arbitrarily  high  accuracy  may 
be  attained  by  allowing  the  step  sizes  to  beccme  small  enough.  While  it  neglects  rounding  error 
in  the  computer  which  becomes  significant  as  many  steps  are  taken,  whether  this  feature  holds 
is  a gauge  on  the  value  of  a difference  scheme.  The  scheme  is  said  to  be  consistent  with  the 
differential  equation  and  its  boundary  conditions  if  the  solution  to  the  approximating  differ- 
ence equations  converges  to  the  exact  solution  of  the  differential  equation  as  the  step  sizes 
vanish.  The  rate  of  convergence  is  determined  by  the  truncation  error  which  is  the  difference 
between  approximate  solution  and  the  exact  solution. 

The  exact  solution  u and  its  coefficients  d!^'  satisfy  (24)  with  M = “>,  N = »>  and 

13  iin 

X = 1/2,  whereas  the  computed  solution  u^^  and  its  coefficients  satisfy  (24)  for  finite 

(0) 

M and  N and  \ = 1/2.  The  necessarily  finite  since  the  heat  equation  (1)  is 

known  to  have  analytic  solutions.  We  let  the  truncation  errors  u. . and  d. . be  defined  by 

u.  . = u.  . - uf?’  , d.  . = d.  . - df?’  . 

13  ID  ID  iD^  iDn  iDn 


Subtraction  of  the  two  forms  of  (24)  gives  the  following  set  of  equations  that  governs  the 
errors. 


(30a) 


N d ® d*®^ 

r il,m+2n  , , ,n  _ r il ,m+2n  , 

^ n!  ~ ^ n!  ■ 

n=0  n=N+l 


/ i = l,2,-.-,L 

k)"  I 

vm  = 0,1, ,M-1 


(30b) 


.(e) 


/'  M <3-  . 00  'N 

^ Jj^2n  ^_^jm  ^ ^ (_j,)m 

^ m ! ^ m ! 


m=0 


m=M+l 


y ,m-K2n  ^ _ y “Li  ,m+2n  j^m 
(•ml  (•  m! 


V m=0 


m=M+l 


j = 1,2,  • • • ,“ 
n ® 0 , 1 , • • • ,N-1 


(30c) 


N d.  . (-k)  - d.  . , , k 

I3,m+2n  i,3-l,m-t-2n 

n=0 


n! 


.(e)  , , ,n  ,(e) 

® d. . - (-k)  - d.  , , 

13  ,m-f2n  i,3-l,m+2n 


n fi  = 1,2,--*,L 


I 

n*N+l 


n! 


M d. . ^ (-h) 

^ J-3  >m+n 


d . h 

1-1 » j ,in+n 


d!®’  (-h)' 

1^ ,m+n 


. h” 
]■-!,]  ,nn-n 


i = 2, 3,  • ■ • ,L 
j = 1,2,  • • • ,“ 
n = 0,1 , • • • ,2N-1 


M N d,  . ^ + d;  , 

' , . , V r in,ra+2n  n ,m+2n  m,  n 

u..(x.,t.)  = ) > ^ ; — ; — ^ x.t. 

1]  1 ] m!n!  i i 

m=0  n=0 


(d  + d*®'  ) » “ d^®* 

ij,M+2W  i],M+2N  M N _ r r i1 ,m+2n  m n 

M!N!  j ^ £ m!n!  *i  j 


i = 1, 2 , • • • ,L 
j = 1,2,  •••,<» 

- M+1  N+1 

Equations  (30)  yield  solutions  for  the  errors  that  are  0(h  ,1c  ) if  the  co- 

efficient matrix  of  the  d. . has  a determinant  that  is  not  zero.  If  this  is  satisfied  then 

i3n 

the  error  in  the  solution  u^^  is  0(h^^^,k*^  and  the  difference  scheme  (24)  is  con- 


sistent with  (1) , (2)  and  (3) . Thus  it  remains  to  examine  the  homogeneous  system. 

In  considering  the  homogeneous  system  we  talce  the  errors  at  the  j-1  st  time  step  to  be 
small  and  investigate  the  error  introduced  due  to  advancing  to  the  jth  time  step.  The  growth 
of  errors  already  introduced  is  a matter  of  stability  and  was  discussed  in  section  4.  Deleting 
j,  the  homogeneous  ecruations  are,  for  all  j. 


N d. 

n • 


i = 1,2,'-*  ,L 


m = 0, 1 , • • • ,M-1 


A-t^+Zn  m ^ Q 

m! 


L,m+2n  .m 

: n =0 

m! 


n = 0 # 1 , • • • »N-1 


M d.  (-h)  - d 

r i ,m4-n 

„ m! 

m=0 


■a  n ('i  = 2,3,'*-,L 

1-1, mtn  = 0 ( 

‘ I n = 0,1,*--,2N-1 


The  coefficient  matrix  of  (31)  resembles  the  Sylvester  matrix  (25)  in  that  a large  number 
of  rows  are  repeated  with  shifts  to  the  right.  This  feature  permits  us  to  find  solutions  that 


contain  powers  of  some  number  z that  can  be  factored  out  of  successive  equations.  It  turns 

out  that  to  obtain  all  of  the  indeoendent  solutions  we  must  allow  different  coefficients  for 

even  and  odd  subscripts  n of  d . . Thus  we  set 

i.n 


d.  _ = A.z 

i,2n  1 


<3.  o = B.z 
i,2n+l  i 


n = 0,1, + J - 


n = 0,1,  • • • , [N  + j - 11 


i = 1,2,-  ••  ,L 


The  independent  equations  after  substitution  into  (31)  are 

N 2n 

(33a)  y ^ (-k)"  = 0 

n=0 


'"l^l  - ®l^-2  = ° 


Vl  ^ ®L^2  = ° 


A.Z,  - B.Z_  - A.  ,Z,  - B.  ,Z,  = 0 
1 1 1 2 1-1  1 1-1  2 


A.Z.,  + B.Z,  - A.  ,Z,  - B.  ,Z,  = 0 
1 2 1 1 1-1  2 1-1  1 


i = 2,3, • • • ,L 


^ y izni  _ 

1 (2m)! 

m=0 


2m  2 

— and  Z = I 


2 (2m+l)! 

in=0 


and  and  are  defined  by  (17). 

There  are  non-trivial  solutions  to  (31)  if  there  are  non-trivial  solutions  for  A.  and 


with  z satisfying  both  (33a)  and  (33b).  In  turn,  there  are  non-trivial  solutions  to  (33b) 
if  the  2L  X 2l  matrix  of  coefficients  is  zero  and  this  condition  is 


»IIHIIUIU'|Plll 


0 0 
0 0 


Equations  (33a)  and  (34)  , for  a fixed  choice  of  M,  N and  L,  give  a set  of  curves  of  Y. 
against  h with  parsuneter  z that  must  be  avoided  vrfien  choosing  h and  Jc  in  order  that  the 
errors  in  (30)  are  of  the  orders  stated.  It  would  seem  improbable  that  a random  choice  of  h 
and  )c  would  satisfy  these  equations.  In  addition,  these  equations  appear  to  be  unrelated  to 
the  stability  conditions  of  section  4,  hence  a choice  of  h and  )c  that  gives  a stable  scheme 
has  a good  chance  of  giving  a consistent  scheme . A choice  of  h and  )c  that  always  produces 
a consistent  scheme  is  found  as  follows.  Choose  h,  and  setting  = Z^,  solve  for  z and 
substitute  into  (33a)  to  get  )c.  Then,  while  (33a)  is  satisfied,  (34)  is  not  since  all  the 
rows  of  the  matrix  are  mutually  orthogonal. 

We  have  not  yet  shown  that  substitution  (32)  gives  all  the  solutions  to  (31)  . This  is  a 
matter  of  chec)ting  that  the  number  of  solutions  for  )c  given  h,  determined  by  setting  the 
determinant  of  the  coefficient  matrix  of  (31)  equal  to  zero,  is  the  same  as  the  number  of 
solutions  for  )c  resulting  from  solving  (34)  for  z and  substituting  into  (33a)  . This  is  a 
difficult  tas)t  because  the  degree  of  eouation  (34)  is  not  easily  found  in  the  general  case.  It 
has  been  chec)ced  in  a large  number  of  examples,  however,  and  it  appears  that  (32)  does  indeed 
lead  to  all  solutions  of  (31) . 


6.  Lowest  Order  Scheme  and  Concluding  Remarks 

If  we  set  M = 1,  N = 1,  X = 1/2  in  (24)  and  eliminate  the  d. . in  favor  of  the  u. . 

tin  11 

we  obtain  the  following  difference  scheme,  arranged  to  display  the  approximation  of  the  deriva- 
tives of  (1)  . 


4 (2k) 


It  has  the  Crank-Nicolson  flavor,  the  approximation  of  the  second  derivative  with  respect  to  x 
being  an  average  of  second  derivative  formulas  at  the  j-1  st  and  jth  time  steps  and  the 
approximation  of  the  time  derivative  being  a weighted  average  over  the  i-1  st,  ith  and  i+1  st 
steps.  This  is  in  fact  the  Keller  Box  Scheme  [4)  for  the  heat  equation  and  can  be  found  by 
writing  (1)  as  a first  order  system  and  using  Keller's  formulas  for  the  first  derivatives.  In 
[4]  Keller  discusses,  for  parabolic  equations  that  include  the  heat  equation,  A-stability,  con- 
vergence, Richardson  extrapolation  and  methods  of  solving  the  difference  equations,  some  of 
which  material  will  carry  over  to  the  higher  order  schemes.  We  speculate  that  the  lowest  order 
scheme  in  the  power  series  method  for  X = 1/2,  applied  to  any  partial  differential  equation, 
coincides  with  Keller's  Box  Scheme  for  that  ecniation  with  a mild  modification  for  problems  in- 
volving free  surfaces.  We  intend  to  show  elsewhere  that  the  power  series  method  can  treat 
irregular  domains  and  free  surface  problems  as  a matter  of  course  and  that  the  lowest  order 
scheme  coincides  with  Keller's  Box  Scheme  with  an  additional  equation  that  adjusts  the  back- 
ground grid  of  cells  to  account  for  the  changing  domain  boundaries.  In  effect,  a co-ordinate 
transformation  is  introduced  automatically  by  the  power  series  method  that  i,.  other  finite 
difference  iiiet)iods  must  be  made  before  approximatipg  the  equations. 

In  conclusion,  the  power  series  n.  hod  has  been  applied  to  the  heat  equation  and  a very 
accurate  difference  representation  has  been  produced.  The  procedure  followed  here  suggests 
that  one  should  discretize  in  each  variable  separately  and  accumulate  stability  conditions. 

In  this  example  discretization  in  t indicated  that  the  point  of  expansion  within  each  cell 
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should  be  on  the  vertical  mid-line,  above  the  center  of  the  cell.  Then  any  accuracy  in  t is 
acceptckble  for  mild  restrictions  on  the  half  time  step  k.  Discretization  in  x indicated 
that  the  scheme  is  stable  for  orders  of  accuracy  in  x less  than  16  for  mild  restrictions  on 
the  half  space  step  h.  The  manner  of  truncation  in  each  variable  indicated  how  to  obtain  a 
well  posed  system  when  both  the  series  are  truncated  simultaneously.  The  stability  conditions 
under  total  discretization  were  modified  to  require  that  the  point  of  expansion  be  precisely 
the  center  of  the  cell  for  mild  restrictions  on  h and  k,  while  the  condition  that  M £ 15 
persisted.  The  "mild  restrictions"  just  referred  to  are  geometrical  and  derive  from  the  root 
locus  diagram  of  the  key  equation  (1-16)  that  was  involved  in  proving  A-stability  of  the  method. 
We  shall  show  in  a succeeding  paper  in  this  series  by  an  analogous  treatment  of  the  wave  equa- 
tion that  the  stability  conditions  of  the  semi-discrete  schemes  may  actually  be  more  restrictive 
than  for  the  totally  discrete  scheme.  It  turns  out  that  one  can  obtain  an  arbitrarily  accurate 
totally  discrete  scheme  for  the  wave  equation  but  not  arbitrarily  accurate  semi-discrete 
schemes  of  the  power  series  type.  This  treatment  also  depends  heavily  on  equation  (1-16) . 
Finally,  the  condition  that  the  scheme  be  consistent  with  the  partial  differential  equation  and 
the  boundary  conditions  provides  further  restrictions  on  h and  k,  and  when  the  scheme  is 
consistent  the  truncation  errors  decrease  as  the  orders  of  accuracy  in  x and  t increase. 
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